Computational complexity and fundamental limitations to fermionic quantum Monte 

Carlo simulations 
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Quantum Monte Carlo simulations, while being efficient for bosons, suffer from the "negative sign 
problem" when applied to fermions - causing an exponential increase of the computing time with 
the number of particles. A polynomial time solution to the sign problem is highly desired since it 
would provide an unbiased and numerically exact method to simulate correlated quantum systems. 
Here we show, that such a solution is almost certainly unattainable by proving that the sign problem 
is NP-hard, implying that a generic solution of the sign problem would also solve all problems in 
the complexity class NP (nondeterministic polynomial) in polynomial time. 



43 

o 



CO 



I 

C 

O 

o 



> 
o 

m 
oo 
o 

o 



-a 
c 

o 
o 



X 



Half a century after the seminal paper of Metropolis 
et al. 0] the Monte Carlo method has widely been estab- 
lished as one of the most important numerical methods 
and as a key to the simulation of many-body problems. 
Its main advantage is that it allows phase space integrals 
for many-particle problems, such as thermal averages, to 
be evaluated in a time that scales only polynomially with 
the particle number N although the configuration space 
grows exponentially with N. This enables the accurate 
simulation of large systems with millions of particles. 

Monte Carlo simulations of quantum systems, such as 
fermions, bosons, or quantum spins, can be performed 
after mapping the quantum system to an equivalent clas- 
sical system. For fermionic or frustrated models this 
mapping may yield configurations with negative Boltz- 
mann weights, resulting in an exponential growth of 
the statistical error and hence the simulation time with 
the number of particles, defeating the advantage of the 
Monte Carlo method. A polynomial time solution of this 
"sign problem" of negative weights would revolutionize 
electronic structure calculations by providing an unbi- 
ased and approximation-free method to study correlated 
fermionic systems. This would be of invaluable help, for 
example, in finding the mechanism for high-temperature 
superconductivity or in determining the properties of 
dense nuclear matter and quark matter. 

The difficulties in finding polynomial time solutions to 
the sign problem are reminiscent of the apparent impos- 
sibility to find polynomial time algorithms for nondeter- 
ministic polynomial (NP)-complete decision problems, 
which could be solved in polynomial time on a hypothet- 
ical non-deterministic machine, but for which no polyno- 
mial time algorithm is known for deterministic classical 
computers. A hypothetical non-deterministic machine 
can always follow both branches of an if-statement si- 
multaneously, but can never merge the branches again. 
It can, equivalently, be viewed as having exponentially 
many processors, but without any communication be- 
tween them. In addition, it must be possible to check 
a positive answer to a problem in NP on a classical com- 
puter in polynomial time. 



Many important computational problems in the com- 
plexity class NP, including the traveling salesman prob- 
lem and the problem of finding ground states of spin 
glasses have the additional property of being NP-hard, 
forming the subset of NP-complete problems, the hard- 
est problems in NP. A problem is called NP-hard if 
any problem in NP can be mapped onto it with poly- 
nomial complexity. Solving an NP-hard problem is thus 
equivalent to solving any problem in NP, and finding a 
polynomial time solution to any of them would have im- 
portant consequences for all of computing as well as the 
security of classical encryption schemes. In that case all 
problems in NP could be solved in polynomial time, and 
hence NP=P. 

As no polynomial solution to any of the NP-complete 
problems was found despite decades of intensive research, 
it is generally believed that NP^P and no deterministic 
polynomial time algorithm exists for these problems. The 
proof of this conjecture remains as one of the unsolved 
millennium problems of mathematics for which the Clay 
Mathematics Institute has offered a prize of one million 
US$ 3]. In this Letter we will show that the sign problem 
is NP-hard, implying that unless the NP^P conjecture 
is disproven there exists no generic solution of the sign 
problem. 

Before presenting the details of our proof, we will give a 
short introduction to classical and quantum Monte Carlo 
simulations and the origin of the sign problem. In the 
calculation of the phase space average of a quantity A, 
instead of directly evaluating the sum 



(A) = ±$Xc)p(c), Z = J>( C ), 
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over a high-dimensional space f2 of configurations c, a 
classical Monte Carlo method chooses a set of M config- 
urations {ci} from f2, according to the distribution p(ci). 
The average is then approximated by the sample mean 
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within a statistical error 

where Var A is the variance of A and the integrated auto- 
correlation time ta is a measure of the autocorrelations 
of the sequence {A(ci)}. In typical statistical physics ap- 
plications, p(c) — exp(— /3-E(c)) is the Boltzmann weight, 
[3 = 1/fcgT is the inverse temperature, and E(c) is the 
energy of the configuration c. 

Since the dimension of configuration space Q grows lin- 
early with the number N of particles, the computational 
effort for the direct integration Eq. scales exponen- 
tially with the particle number N. Using the Monte Carlo 
approach the same average can be estimated to any de- 
sired accuracy in polynomial time, as long as the auto- 
correlation time ta does not increase faster than polyno- 
mially with N. 

In a quantum system with Hamilton operator H, in- 
stead of an integral like Eq. JQ), an operator expression 



(A) = -Tr[Aexp(-/3#)] - Z 



Trexp(-/?£0 (3) 



needs to be evaluated in order to calculate the thermal 
average of the observable A (represented by a self-adjoint 
operator). Monte Carlo techniques can again be applied 
to reduce the exponential scaling of the problem, but only 
after mapping the quantum model to a classical one. One 
approach to this mapping 4] is a Taylor expansion [f| : 



Trexp(-/iff) = ^^TxH 11 
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where for each order n in the expansion we insert n sums 
over complete sets of basis states {\i)}. The "configura- 
tions" are sequences c = (ix,...,i n ) of n basis states and 
we define the weight p(c) by the corresponding product 
of matrix elements of H and the term (— (3) n /n\. With a 
similar expansion for Tr[Aexp(— (3H)\ we obtain an ex- 
pression reminiscent of classical problems: 
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|Tr[Aexp(-/3^)] = i^^(c)p(c) 
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If all the weights p(c) are positive, standard Monte 
Carlo methods can be applied, as it is the case for non- 
frustrated quantum magnets and bosonic systems. In 
fermionic systems |(| negative weights p(c) < arise 
from the Pauli exclusion principle, when along the se- 
quence — > |«2) \i n ) — * two fcrmions are 
exchanged, as shown in Fig. 1. 

The standard way of dealing with the negative weights 
of the fermionic system is to sample with respect to 




FIG. 1: A configuration of a fermionic lattice model on a 4- 
site square. The configuration has negative weight, since two 
fermions are exchanged in the sequence — > \i<z) — » — » 
— * World lines connecting particles on neighboring 
slices are drawn as thick lines. 



the bosonic system by using the absolute values of the 
weights \p(c)\ and to assign the sign s(c) = signp(c) to 
the quantity being sampled: 



(A) = 



E c A(c) S (c)\p(c)\/j:Mc)\ 

E "Wb(c)l/Ec b( c )l 



(6) 



(As)' 



While this allows Monte Carlo simulations to be per- 
formed, the errors increase exponentially with the par- 
ticle number N and the inverse temperature (3. To see 
this, consider the mean value of the sign (s) = Z/Z 1 , 
which is just the ratio of the partition functions of the 
fermionic system Z = ^2 c p(c) with weights p(c) and the 
bosonic system used for sampling with Z' = ^2 c \p(c)\. 
As the partition functions are exponentials of the cor- 
responding free energies, this ratio is an exponential of 
the differences A/ in the free energy densities: (s) = 
Z/Z' = exp(—f3NAf). As a consequence, the relative 
error As/(s) increases exponentially with increasing par- 
ticle number and inverse temperature: 
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Similarly the error for the numerator in Eq. (7) in- 
creases exponentially and the time needed to achieve a 
given relative error scales exponentially in N and /3. 

In order to avoid any misconception about what would 
constitute a "solution" of the sign problem, we start by 
giving a precise definition: 

• A quantum Monte Carlo simulation to calculate a 
thermal average (A) of an observable A in a quan- 
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turn system with Hamilton operator H is denned 
to suffer from a sign problem if there occur negative 
weights p(c) < in the classical representation as 
given by Eq. JSJ. 

• The related bosonic system of a fermionic quantum 
system is defined as the system where the weights 
p(c) are replaced by their absolute values |p(c)|, 
thus ignoring the minus sign coming from fermion 
exchanges: 

(AY = ^J2 A ( c Mc)\- (8) 

c 

• An algorithm for the stochastic evaluation of a ther- 
mal average such as Eqns. (0 or © is defined to 
be of polynomial complexity if the computational 
time t(e,N,{3) needed to achieve a relative statis- 
tical error e — A A/ (A) in the evaluation of the 
average (A) scales polynomially with the system 
size N and inverse temperature (3, i.e. if there exist 
integers n and m and a constant k < oo such that 

t(e,N,P) < Ke- 2 N n f3 m . (9) 

• For a quantum system that suffers from a sign prob- 
lem for an observable A, and for which there exists 
a polynomial complexity algorithm for the related 
bosonic system Eq. (|SJ| , we define a solution of the 
sign problem as an algorithm of polynomial com- 
plexity to evaluate the thermal average (A) . 

It is important to note that we only worry about the 
sign problem if the bosonic problem is easy (of polyno- 
mial complexity) but the fermionic problem hard (of ex- 
ponential complexity) due to the sign problem. If the 
bosonic problem is already hard, e.g. for spin glasses 0, 
the sign problem will not increase the complexity of the 
problem. Also, changing the representation so that the 
sum in Eq. (jSJ contains only positive terms p(c) > is not 
sufficient to solve the sign problem if the scaling remains 
exponential, since then we just map the sign problem 
to another exponentially hard problem. Only a polyno- 
mial complexity algorithm counts as a solution of the sign 
problem. 

At first sight such a solution seems feasible since the 
sign problem is not an intrinsic property of the quan- 
tum model studied but is representation-dependent: it 
depends on the choice of basis sets and in some 

models it can be solved by a simple local basis change || . 
Indeed, when using the eigenbasis in which the Hamilton 
operator H is diagonal, there will be no sign problem. 
This diagonalization of the Hamilton operator is, how- 
ever, no solution of the sign problem since its complexity 
is exponential in the number of particles N . 

We now construct a quantum mechanical system for 
which the calculation of a thermal average provides the 



solution for one and thus all of the NP-complete prob- 
lems. This system exhibits a sign problem, but the re- 
lated bosonic problem is easy to solve. Since, for this 
model, a solution of the sign problem would provide us 
with a polynomial time algorithm for an NP-complctc 
problem, the sign problem is NP-hard. Of course, it is 
expected that the corresponding thermal averages cannot 
be calculated in polynomial time and the sign problem 
thus cannot be solved. Otherwise we would have found 
a polynomial time algorithm for the NP-complete prob- 
lems and would have shown that NP=P. 

The specific NP-complete problem we consider Q is to 
determine whether a state with energy less than or equal 
to a bound Eq exists for a classical three-dimensional 
Ising spin glass with Hamilton function 

H = -^J jk a j a k . (10) 
U-k) 

Here the spins <jj take the values ±1, and the couplings 
Jjk between nearest neighbor lattice points j and k are 
either or ± J. 

This problem is in the complexity class NP since the 
non-deterministic machine can evaluate the energies of 
all configurations c in polynomial time and test whether 
there is one with E{c) < Eq. In addition, the validity 
of a positive answer (i.e. there is a configuration c) can 
be tested on a deterministic machine by evaluating the 
energy of that configuration. The evaluation of the par- 
tition function Z = ^ c exp(— (3E{c)) is, however, not in 
NP since the non-deterministic machine cannot perform 
the sum in polynomial time. 

This question whether there is a state with energy 
E(c) < Eq can also be answered in a Monte Carlo sim- 
ulation by calculating the average energy of the spin 
glass at a large enough inverse temperature f3. Since 
the energy levels are discrete with spacing J it can eas- 
ily be shown that by choosing an inverse temperature 
(3 J > N In 2 + ln(12iV) the thermal average of the energy 
will be less than Eq + J/2 if at least one configuration 
with energy Eq or less exists, and larger than Eq + J 
otherwise |9J. 

In this classical Monte Carlo simulation, the complex 
energy landscape, created by the frustration in the spin 
glass (Fig. 2a), exponentially suppresses the tunneling 
of the Monte Carlo simulation between local minima at 
low temperatures. The autocorrelation times and hence 
the time complexity of this Monte Carlo approach are 
exponentially large t oc exp(aA), as expected for this 
NP-complete problem. 

We now map this classical system to a quantum system 
with a sign problem. We do so by replacing the classical 
Ising spins by quantum spins. Instead of the common 
choice in which the classical spin configurations are basis 
states and the spins are represented by diagonal cr| Pauli 
matrices we choose a representation in which the spins 
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FIG. 2: a) A classically frustrated spin configuration of three 
antiferromagnetically coupled spins: no configuration can si- 
multaneously minimize the energy of all three bonds, b) A 
configuration of a frustrated quantum magnet with negative 
weights: three antiferromagnetic exchange terms with nega- 
tive weights are present in the sequence —> |i2) — > |i3) — ► 
Here up-spins with ^-component of spin cr| = 1 and 
down-spins with cr| = — 1 are connected with differently col- 
ored world lines. 

point in the ±x direction and are represented by Pauli 
matrices: 

H=-J2jjk<T]<r%, (11) 

Here the random signs of the couplings are mapped to 
random signs of the off-diagonal matrix elements which 
cause a sign problem (see Fig. 2b). The related bosonic 
model is the ferromagnet with all couplings Jjk > and 
efficient cluster algorithms with polynomial time com- 
plexity are known for this model 10] . Since the bosonic 
version is easy to simulate, the sign problem is the origin 
of the NP -hardness of a quantum Monte Carlo simula- 
tion of this model. A generic solution of the sign problem 
would provide a polynomial time solution to this, and 
thus to all, NP-complete problems, and would hence 
imply that NP=P. Since it is generally believed that 
NP^P, we expect that such a solution does not exist. 

Conclusions - By constructing a concrete model we 
have shown that the sign problem of quantum Monte 
Carlo simulations is NP-hard. This does not exclude 
that a specific sign problem can be solved for a restricted 
subclass of quantum systems. This was indeed possible 
using the meron-cluster algorithm for some partic- 
ular lattice models. Such a solution must be intimately 
tied to properties of the physical system and allow an 
essentially bosonic description of the quantum problem. 
A generic approach like the cancellation idea 12] might 
scale polynomially for some cases but will in general scale 
exponentially. 

In the case of fermions or frustrated quantum mag- 
nets, solving the sign problem requires a mapping to a 



bosonic or non-frustrated system - which is, in general, 
almost certainly impossible for physical reasons. The ori- 
gin of the sign problem is, in fact, the distinction between 
bosonic and fermionic systems. The brute-force approach 
of taking the absolute values of the probabilities means 
trying to sample a frustrated or fermionic system by sim- 
ulating a non-frustrated or bosonic one. As for large 
system sizes N and low temperatures the relevant con- 
figurations for the latter are not the relevant ones for the 
former, the errors are exponentially large. 

Given the NP-hardness of the sign problem one 
promising idea for the simulation of fermionic systems 
is to use ultra-cold atoms in optical lattices to construct 
well-controlled and tunable implementations of physical 
systems, such as the Hubbard model 0, and to use 
these "quantum simulators" to study the phase diagrams 
of correlated quantum systems. But even these quantum 
simulators are most likely not a generic solution to the 
sign problem since there exist quantum systems with ex- 
ponentially diverging time scales and it is at present not 
clear whether a quantum computer could solve the NP- 
complete problems [lij . 
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